###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*##
###*### COMPUTE COUNTERFACTUAL EXPERIMENTS              ###*###*###*###*###*###*###*###*##
###*### USING PREVIOUSLY ESTIMATED PARAMETERS           ###*###*###*###*###*###*###*###*##
###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*##

#' FROM THE PAPER:
#' Child Care Markets, Parental Labor Supply, and Child Development
#' BY:
#' Samuel Berlinski, Maria Marta Ferreyra, Luca Flabbi and Juan David Martin (2023)

# ___________________________________________________________________________________ ####
# START ####

## External packages
library(lubridate)
library(snowfall)
library(readxl)

rm(list = ls())
gc()

## Directory
folderpath <- '~/EstimationCode/'
setwd(folderpath)

# ___________________________________________________________________________________ ####
# INPUTS ####

## Functions
source('functions.R')
source('simul_market.R')
source('simul_market_PP.R')

## Simulated ex ante data
load('SimData_M0.RData')

## Parameter estimates
load('Estimates_01.RData')

# ___________________________________________________________________________________ ####
# SETUP ####

## Market given settings
J     <- nrow(data)           # number of hh types
NL    <- 7                    # No. of potential entrants
M     <- 40*NL/J/0.3          # Total market size
Fixed <- c(128.90, 322.25)    # Fixed cost for LQ and HQ, respectively
minw  <- 5.15                 # Minumum wage
Cn    <- minw                 # Opp. cost of non-relative caregiver
Cr    <- 0                    # cost of relative caregiver
I     <- 41                   # daily non-labor income
subs  <- 0.35                 # subsidy discount for single moms
Nrel  <- round(J*0.4)         # No. of hh with relative care available

## Observed heterogeneity
q0  <- data$q0      # mental score age 1
ed  <- data$edu     # mom's edu
edf <- data$edf     # dad's edu
wf  <- data$wf      # dad's wage
sm  <- data$single  # single moms
To  <- 16 + 2*sm    # time endowment

## Dad's wage implied random term
ef  <- log(wf)      
mwf <- mean(ef[sm == 0])
swf <- sd(ef[sm == 0])
ef  <- (ef - mwf)/swf
ef[wf == 0] <- 0

## Random draws
set.seed(123)
ev   <- matrix(rnorm(3*J), J, 3) # preferences
ew   <- matrix(rnorm(J),   J, 1) # wages
ec   <- rowMeans(apply(matrix(rnorm(NL*10000), NL, 10000), 2, sort)) # marginal costs
ec   <- (ec - mean(ec))/sd(ec)
reli <- sample(1:J, size = Nrel, replace = F) # hh with relative care available

## Standarization for wages draws
ew[ed == 1] <- (ew[ed == 1] - mean(ew[ed == 1]))/sd(ew[ed == 1])
ew[ed == 2] <- (ew[ed == 2] - mean(ew[ed == 2]))/sd(ew[ed == 2])
ew[ed == 3] <- (ew[ed == 3] - mean(ew[ed == 3]))/sd(ew[ed == 3])

# ___________________________________________________________________________________ ####
# SIMULATE BASELINE ####

## Simulate
keep <- c(objects(), "keep")
sfInit(parallel = TRUE, cpus = 4)
sfExport(list = keep, local = FALSE)
market <- simul_market(pp1)
sfStop()

## Save
market_BL <- market
save(market_BL, file = 'simBL.RData')

# ___________________________________________________________________________________ ####
# SIMULATE CF1: CASH TRANSFER ####

## Simulate
keep <- c(objects(), "keep")
sfInit(parallel = TRUE, cpus = 4)
sfExport(list = keep, local = FALSE)
market <- simul_market(pp1, Cash = 40)
sfStop()

## Save
market_CF1 <- market
save(market_CF1, file = 'simCF1.RData')

# ___________________________________________________________________________________ ####
# SIMULATE CF2: VOUCHER UNCOND. ####

## Simulate
keep <- c(objects(), "keep")
sfInit(parallel = TRUE, cpus = 4)
sfExport(list = keep, local = FALSE)
market <- simul_market(pp1, tau = 5)
sfStop()

## Save
market_CF2 <- market
save(market_CF2, file = 'simCF2.RData')

# ___________________________________________________________________________________ ####
# SIMULATE CF3: VOUCHER MOM WORK ####

## Simulate
keep <- c(objects(), "keep")
sfInit(parallel = TRUE, cpus = 4)
sfExport(list = keep, local = FALSE)
market <- simul_market(pp1, tau = 5, both = T)
sfStop()

## Save
market_CF3 <- market
save(market_CF3, file = 'simCF3.RData')

# ___________________________________________________________________________________ ####
# SIMULATE CF4: VOUCHER INCOME ####

## Simulate
keep <- c(objects(), "keep")
sfInit(parallel = TRUE, cpus = 4)
sfExport(list = keep, local = FALSE)
market <- simul_market(pp1, tau = 5, both = F, mini = 90)
sfStop()

## Save
market_CF4 <- market
save(market_CF4, file = 'simCF4.RData')

# ___________________________________________________________________________________ ####
# SIMULATE CF5: VOUCHER QUALITY ####

## Simulate
keep <- c(objects(), "keep")
sfInit(parallel = TRUE, cpus = 4)
sfExport(list = keep, local = FALSE)
market <- simul_market(pp1, tau = 5, hqv = T)
sfStop()

## Save
market_CF5 <- market
save(market_CF5, file = 'simCF5.RData')

# ___________________________________________________________________________________ ####
# SIMULATE CF6: QUALITY REQ. ####

## Simulate
keep <- c(objects(), "keep")
sfInit(parallel = TRUE, cpus = 4)
sfExport(list = keep, local = FALSE)
market <- simul_market(pp1, hqr = T)
sfStop()

## Save
market_CF6 <- market
save(market_CF6, file = 'simCF6.RData')

# ___________________________________________________________________________________ ####
# SIMULATE CF7: PUBLIC PROVISION P = MeC ####

## Simulate
keep <- c(objects(), "keep")
sfInit(parallel = TRUE, cpus = 4)
sfExport(list = keep, local = FALSE)
market <- simul_market_PP(pp1)
sfStop()

## Save
market_CF7 <- market
save(market_CF7, file = 'simCF7.RData')

# ___________________________________________________________________________________ ####
# END ####